###################################
###### This script contains examples of how to fit distributions to datasets

# Load library

library('SDMTools')

## Identify directory for ASCII files

asc.dir = '/home1/99/jc152199/MAXENT/output/'

# Identify species of interest

species = c('SAPCZEC','SAPTETR','CARRUBR','LAMCOGG','GNYQUEE','SAPBASI','LAMROBE')

## Blank object for data binding within loop

outdata = NULL

### Loop through species, extracting and summarizing ClassStat data

for (spp in species)

	{
	
	### Loop through different models within species, extracting and summarizing data
	
	for (mod in c('BIOCLIM','microCLIM'))
	
		{
		
		### Read in the ASCII
		
		t.asc = read.asc(paste(asc.dir,spp,'/',mod,'/output/',spp,'_',mod,'.asc',sep=''))
		
		### Remove NA's
		
		t.asc = t.asc[which(is.finite(t.asc)==T)]
		
		### Creating a histogram
		
		hist(t.asc,main="Hist of All ES Values")
		
		dev.off()
		
		######## Plot density
		
		plot(density(t.asc),main="Density estimate of data")
		
		dev.off()
		
		#### Plot empirical cumulative distribution function
		
		plot(ecdf(t.asc),main= 'Empirical cumulative distribution function')
		
		dev.off()
		
		#### Q-Q Plot (based on normally distributed data I think)
		
		png('check.png')
		z =(t.asc-mean(t.asc))/sd(t.asc) ## standardized data
		qqnorm(z) ## drawing the QQplot
		abline(0,1) ## drawing a 45-degree reference line
	
		dev.off()
		
		### Fitting to a Weibull distribution
		
		
		x=seq(0,1,.1)
		dgamma(x, scale=1.5, shape=2)
		curve(dgamma(x, scale=1.5, shape=2),from=0, to=1)
		
		}
		
	}
	
### Done